function eq_mit = solve_eq_mit(sim, eq_SS, param, glob, options)

    diff = 10;
    %iter = 0;
    
    T_plot = 30;
    
    
    ct = 0;
    
    while (diff > 1e-6)
    
        ct = ct+1;
        
        sim_out  = solve_bf(sim, eq_SS, param, glob, options);
        
        diffC = max(max(abs(sim_out.C - sim.C)/max(abs(sim.C))));
        diffD = max(max(abs(sim_out.D - sim.D)/max(abs(sim.D))));
        diffW = max(max(abs(sim_out.W - sim.W)/max(abs(sim.W))));
        
        diff  = max([diffC diffD diffW])
        
%         Labor       =  sim_out.l;
%         total_labor =  sim_out.L; sum(sim_out.mu.*Labor);
%         
%         figure(47)
%         
%         subplot(3,3,1)
%         hold off
%         plot([sim_out.C - 1]); drawnow; title('C'); drawnow;
%         xlim([1 T_plot]);
%         
%         figure(47)
% 
%         subplot(3,3,2)
%         hold off
%         plot(sim_out.D - eq_SS.D); drawnow; title('D'); drawnow;
%         xlim([1 T_plot]);
% 
%         figure(47)
%         
%         subplot(3,3,3)
%         hold off
%         plot([total_labor(end), total_labor]); title('Total labor demand'); drawnow;
%         xlim([1 T_plot]);
% 
%         figure(47)
%         
%         subplot(3,3,4)
%         hold off
%         plot(sim_out.entry_mass./sum(sim_out.mu)); title('Entry rate'); drawnow;
%         xlim([1 T_plot]);
%         
%         figure(47)
%         subplot(3,3,5)
%         hold off
%         plot(sum(sim_out.exit.*sim_out.mu)./sum(sim_out.mu)); title('Exit rate'); drawnow;
%         xlim([1 T_plot]);
% 
%         figure(47)
%         subplot(3,3,6)
%         hold off
%         plot(sim_out.employment_young./sim_out.L); title('Share of L at firms < 5yo'); drawnow;
%         xlim([1 T_plot]);
% 
%         costs = sim_out.l;
%         prods = repmat(exp(glob.sf(:,2)), 1, options.T);
%         wages = repmat(sim_out.W', glob.Nsf, 1);
%         mkps  = sim_out.p./(wages./prods);
%         
%         tc    = sum(costs.*sim_out.mu);
%         num    = sum(costs.*sim_out.mu.*mkps);
%         
%         mu_cw = num./tc;
%         
%         figure(47)
%         subplot(3,3,7)
%         hold off
%         plot((mu_cw)); title('Cost-weighted markup'); drawnow;
%         xlim([1 T_plot]);
% 
%         
%         figure(47)
%         subplot(3,3,8)
%         hold off
%         plot((sim_out.W)); title('Wage');  drawnow;
%         xlim([1 T_plot]);
% 
%         % plot share of employment among firms <1 years old
%         figure(47)
%         subplot(3,3,9)
%         hold off
%         plot(100*sim_out.entrant_employment./total_labor);
%         title('Pct share of employment among entrants');
%         xlim([1 10]);
%         drawnow;
%  
%         if diff < .025
%             options.damp_mit = .99;
%         end
%         
%         if diff < 0.01
%             options.damp_mit = .995;
%         end
%         
        figure(47)
        subplot(4,1,1)
        hold off
        plot(sim.C(1:30));
        hold on
        plot(sim_out.C(1:30)); 
        title('C'); legend('in', 'out'); drawnow;
        
        figure(47)
        subplot(4,1,2)
        hold off
        plot(sim.D(1:30));
        hold on
        plot(sim_out.D(1:30));
        title('D'); legend('in', 'out'); drawnow;
        
        figure(47)
        subplot(4,1,3)
        hold off
        plot(sim.W(1:30));
        hold on
        plot(sim_out.W(1:30));
        title('W'); legend('in', 'out'); drawnow;

        figure(47)
        subplot(4,1,4)
        hold off
        plot(sum(sim_out.mu(:,1:30))/sum(sim_out.mu(:,end)));
        title('Mass'); drawnow;

        
        % recalculate labor
         
        diffC = sim_out.C - sim.C;
        diffC = abs(diffC)/max(abs(diffC));
        diffC = diffC;
        diffC = (1-options.damp_mit)*diffC;
        sim.C =(1-diffC).*sim.C + diffC.*sim_out.C;
        
        diffD = sim_out.D - sim.D;
        diffD = abs(diffD)/max(abs(diffD));
        diffD = diffD;
        diffD = (1-options.damp_mit)*diffD;
        sim.D =(1-diffD).*sim.D + diffD.*sim_out.D;

        diffW = sim_out.W - sim.W;
        diffW = abs(diffW)/max(abs(diffW));
        diffW = diffW;
        diffW = (1-options.damp_mit)*diffW;
        sim.W =(1-diffW).*sim.W + diffW.*sim_out.W;
        
        
        
    end
    
    sim_out.entry = [];
    eq_mit = sim_out;
    
end